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INTRODUCTION 

In  the  present  study,  we  shall  consider  the  problem  of  determining  the  ability  of  fiber  rein¬ 
forced  composite  material  structures  to  sustain  handling  and  impact  loading.  Such  an  investigation 
will  represent  an  important  element  in  assessing  the  applicability  of  such  material  systems  for  long 
term  service  utilization.  Accordingly,  the  present  work  is  directed  toward  the  development  of  a 
technique  which  will  predict  the  strains  in  clamped  plates  which  arise  as  a  consequence  of  relatively 
low  speed  (0  to  30  m/s),  hard  object,  transverse  normal  impact. 

Previously  we  have  considered  the  problems  of  the  low  velocity  impact  of  both  beams  and 
simply-supported  orthotropic  plates,  i.e.,  references  (1),  (2)  and  (3).  Thus,  in  reference  (1)  we 
initiated  our  study  of  low  velocity  impact  by  considering  the  analysis  of  a  simply  supported  beam. 

In  this  first  investigation  it  was  shown  that  the  physical  problem  of  interest  was  related  to  the  struc¬ 
tural  response  of  the  system.  Additionally,  it  was  determined  that  the  destruction  of  the  beam  was 
related  to  a  tension  failure  in  the  extreme  fiber  opposite  to  the  impact  point.  In  reference  (2)  one 
of  the  authors  considered  the  effect  of  the  introduction  of  viscoelastic  damping  into  the  solution  of 
a  simply  supported  beam  problem.  Since  the  undamped  theoretical  solution  tends  to  form  an  upper 
bound  for  the  experimental  data,  this  work  led  to  an  improved  correlation  of  the  theoretical  and 
experimental  data.  Finally,  in  reference  (3)  one  of  the  authors  considered  the  extension  of  the  un¬ 
damped  theory  to  the  problem  of  a  simply  supported  orthotropic  plate.  However,  it  is  important 
to  observe  that  this  last  work  was  extended  to  include  damping  by  empirically  correcting  the  solu¬ 
tion  according  to  the  results  of  reference  (2).  Indeed,  we  may  summarize  the  essence  of  all  of  the 
previous  work  by  declaring  that  an  undamped  theoretical  solution  of  a  vibrating  system  suitability 
corrected  for  damping  correlates  the  available  experimental  data,  reference  (4),  to  within  the  limits 
of  experimental  accuracy. 

In  view  of  the  success  of  the  previous  investigations  it  was  decided  to  extend  the  work  to  the 
practically  important  problem  of  a  clamped  plate.  Since  the  theoretical  solution  of  a  clamped  plate 
problem  is  more  difficult  than  a  corresponding  simply  supported  plate  problem  and  the  solution  of 
the  simply  supported  plate,  i.e.,  reference  (3),  extended  over  three  hundred  pages,  a  new  solution 
technique  was  indicated.  Thus,  we  considered  the  solution  of  a  clamped  plate  subjected  to  impact 
utilizing  the  finite  element  approach.  Although  the  present  problem  resides  within  the  capacities  of 
most  general  purpose  finite  element  programs,  i.e.,  NASTRAN,  it  was  decided  that  a  stand  alone 
program  which  was  directed  at  the  solution  of  this  one  problem  would  allow  a  more  economical 
parametric  study  of  the  general  situation.  Indeed,  an  initial  investigation  of  two  specific  clamped 
plate  problems  employing  NASTRAN  required  an  average  of  thirty-three  seconds  of  computer  CPU 
time  on  the  CDC  computer  which  translated  into  an  approximately  sixty  three  dollar  per  problem 
charge.  Certainly  a  large  parametric  study  subject  to  such  constraints  would  be  prohibitive. 
NASTRAN  permits  the  utilization  of  the  inverse  power  method  in  the  solution  of  the  eigenvalue 
problem,  however,  it  provides  for  no  initial  estimate  of  the  magnitude  of  the  eigenvalue.  Indeed, 
NASTRAN  avoids  this  problem  by  employing  an  iterative  search  technique  which  employs  range 
estimates  supplied  by  the  user.  Thus,  a  stand  alone  program  which  would  supply  initial  estimates 
of  the  eigenvalues  represented  a  desired  element  of  the  overall  solution.  This  element  of  the  total 
task  was  accomplished  by  utilizing  the  energy  method  to  develop  a  technique  which  supplies 
approximations  for  the  first  two  eigenvalues.  Experimental  data  indicated  that  damping  attenuated 
the  contribution  of  the  other  mode  shapes  sufficiently  that  they  were  not  visible  in  the  output 
trace  of  the  experimental  records  in  the  time  domain  under  consideration.  The  program  which 
resulted  from  this  analysis  produced  estimates  of  the  first  and  second  eigenvalues  which  were,  on 
the  average,  ten  and  thirty  percent  larger  than  the  corresponding  finite  element  values.  Since  the 
energy  method  always  yields  approximate  values  for  the  eigenvalues  which  are  greater  than  the 
minimum  possible  value  in  the  space  under  consideration,  the  above  factors  were  employed  to 
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correct  the  initial  frequency  estimates.  Specifically,  the  first  and  second  approximate  frequencies 
were  multiplied  by  factors  of  0.92  and  0.76  before  they  were  input  into  the  finite  element  program. 
The  finite  element  program  was  designed  for  computational  efficiency  and  minimum  storage  re¬ 
quirements.  Indeed,  the  program  permits  a  ten  by  ten  subdivision  of  the  plate,  i.e.,  100  nodes, 
which  corresponds  to  approximately  a  three  hundred  degrees-of-freedom  finite  element  mode.  All 
computational  effort  is  performed  in  core  with  a  result  that  for  the  maximum  permissible  sub 
division  of  the  plate  the  program  requires  only  slightly  over  four  seconds  of  CPU  on  the  CDC  com¬ 
puter.  This  translates  into  a  cost  of  approximately  three  and  one  half  dollars  per  problem  since 
three  problems  were  normally  solved  at  one  time.  Finally,  a  third  computer  program  was  written 
which  would  accept  the  results  of  the  finite  element  solution  and  determine  the  maximum  strain 
at  the  center  of  the  plate.  This  program  permitted  the  amplitudes  of  the  mode  shapes  to  be  cor 
rected  for  damping  in  accord  with  results  of  the  beam  damping  analysis,  reference  (2),  i.e.,  we 
employed  the  procedure  utilized  in  reference  (3). 

This  report  extends  the  analysis  of  composite  structural  impact  damage  by  considering  a 
boundary  value  problem  of  considerable  practical  utility.  It  is  to  be  understood  that  the  finite 
element  solution  coincides  with  the  following  physical  problem.  A  mass  of  small  dimensions, 
possessing  unit  velocity,  is  directed  at  a  plate  which  is  initially  at  rest.  During  impact  the  impactor 
transfers  its  momentum  to  the  region  of  the  plate  with  which  it  is  in  contact.  This  information  is 
utilized  to  ascertain  initial  conditions  for  the  vibrating  system.  The  plate  and  impactor  remain  in 
contact,  over  a  small  area,  during  the  remainder  of  the  first  quarter  of  the  first  cycle  of  the  motion. 
Finally,  the  finite  element  solution  is  empirically  corrected,  on  the  basis  of  viscoelastic  damped 
beam  theory,  to  account  for  the  effects  of  damping.  Thus,  the  solution  presented  in  this  report 
considers  central  impact,  provides  the  finite  element  approximation  for  the  eigenvalues  and  eigen 
functions  of  a  clamped  plate  carrying  a  concentrated  impactor  mass  and  it  is  empirically  altered 
to  account  for  material  damping. 


ANALYSIS 

The  work  performed  during  this  investigation  may  conveniently  be  subdivided  into  two  dis¬ 
tinct  subsections.  We  shall  initially  describe  the  basis  of  the  finite  element  model,  then,  having  sum 
marized  the  method  of  securing  the  eigenvalues  and  eigenvectors  of  the  problem,  proceed  with  the 
analysis  of  the  impact  event. 

For  conceptual  clarity,  we  shall  divide  our  description  of  the  finite  element  model  into  two 
unique  subtasks.  First,  we  shall  describe  the  finite  element  model  and  our  solution  of  the  resulting 
system  of  equations,  then,  we  shall  describe  our  method  of  approximating  the  eigenvalues  of  a 
vibrating  clamped  plate. 

FINITE  ELEMENT  MODEL 


The  potential  energy  of  a  vibrating  plate  may  be  shown  to  be  expressible  in  the  form. 


/ 


U(w)  dA 


where: 


(1) 


A  =  area  of  plate 
PE  =  potential  energy 
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U<w)  -D1lW2x  +  2H  w2y  +  D22  wyy 
D1 1'  d22  •  •  •  =  flexural  rigidities  of  the  plate 
H  =  °12  +  2d66 

w  =  transverse  displacement  of  the  plate 


Note  that  subscripts  appended  to  the  transverse  displacement  indicate  partial  derivatives  with  re¬ 
spect  to  the  appropriate  space  variable  and  a  dot  denotes  a  time  derivative.  Further,  the  kinetic 
energy  of  the  plate  may  be  written  as, 

KE=1  /phw2dA  (2) 

*  A 


where: 


KE  =  kinetic  energy 
h  =  thickness  of  plate 
p  =  mass  density  of  plate 

Since  we  are  concerned  with  the  determination  of  the  eigenvectors,  i.e.,  mode  shapes,  of  the 
vibrating  system  we  will  assume  that  the  motion  is  harmonic,  or, 

w  =  w  sin  cjt  (3) 


where: 


w  =  mode  shape 
t  =  time 


«  =  circular  frequency 


Introducing  the  expression  for  the  mode  shape,  equation  (3),  into  the  kinetic  energy  formula, 
equation  (2),  we  find  that  the  maximum  kinetic  energy,  KEmax,  is  given  by, 

KEmax  =  ,Jw2Ph  /w2dA 
A 


(4) 


The  results  expressed  by  equations  (1)  and  (4)  were  utilized  to  construct  a  finite  element  represen¬ 
tation  of  our  problem.  Details  of  this  work  are  summarized  in  appendix  A. 


APPROXIMATE  NATURAL  FREQUENCIES 


Next  we  turn  our  attention  to  the  problem  of  the  approximate  determination  of  the  first  two 
natural  frequencies  of  the  system.  This  task  is  an  important  element  of  our  procedure  since  the 
determination  of  reasonably  accurate  approximate  frequencies  will  considerably  reduce  computa¬ 
tional  costs.  Indeed,  a  major  portion  of  our  success  in  decreasing  computational  effort  may  be 
traced  to  this  portion  of  our  analysis. 
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In  the  problem  under  consideration  we  recall  that  the  concentrated  mass  impacts  the  origin  of 
the  coordinate  system  and,  therefore,  the  kinetic  energy  of  our  system  may  be  written, 

KE  =^M  Wq  +  1  /*phvv2dA  (5) 

•'A 

where: 


M  =  mass  of  impacting  object 


wg  =  velocity  of  the  plate  at  the  origin  of  the  coordinate  system 


Introducing  the  expression  for  the  mode  shape,  equation  (3),  into  the  potential  energy  expression, 
equation  (1),  and  the  kinetic  energy  expression,  equation  (5),  and  equating  the  resulting  maximums, 
we  obtain, 


L 


a)2  =  — — j 

MWg  +  p  h 


A  U(w)  dA 

h/> 


w2  dA 


(6) 


In  order  to  simplify  our  computations  we  shall  convert  the  variables  under  the  integral  sign  in 
equation  (6)  to  dimensionless  form,  i.e., 


where: 


a  =  half  length  of  plate 
b  =  half  width  of  plate 


In  view  of  the  symmetry  of  the  problem  about  the  origin  the  integrations  were  confined  to  the 
region  of  the  plate  corresponding  to  positive  x  and  y.  This  alteration  in  the  integration  procedure 
requires  that  we  utilize  one  quarter  of  the  impacting  mass  in  equation  (6). 
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The  result  indicated  in  equation  (6)  was  employed  to  determine  approximate  values  of  the  first 
two  natural  frequencies  of  the  vibrating  plate.  All  the  details  of  this  work  are  summarized  in 
appendix  8. 

INVERSE  POWER  METHOD 

Since,  from  the  previous  section,  we  possess  approximations  for  the  first  two  natural  fre¬ 
quencies  of  the  system  it  is  efficient  to  employ  an  iterative  procedure  to  solve  our  finite  element 
equations.  Specifically,  we  shall  employ  the  computational  algorithm  known  as  the  inverse  power 
method. 

After  assembly  the  system  of  finite  element  equations  is  expressible  in  the  form, 

KX  =  \MX  (7) 


where: 


K  =  stiffness  matrix 
M  =  mass  matrix 
X  =  eigenvector 
X  =  eigenvalue 

We  shall  partition  the  eigenvalue  X  into  two  parts,  i.e., 
X  - Xq  +  AX 

where: 


Xq  =  approximation  to  the  eigenvalue 
AX  =  correction  to  the  approximation 


Introducing  the  above  expression  for  the  eigenvalue  into  the  original  set  of  equations,  equations  (7), 
and  normalizing  the  input  vector  so  that  the  largest  component  assumes  the  value  one,  we  obtain, 


(K  -  X0M)X  =  AXXmax 


M 


) 


where: 


Xmax  =  maximum  component  of  the  solution  vector 
However,  it  may  be  shown  that  as  the  number  of  iterations  increases. 


AX  X 


max 


1 
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Therefore,  we  may  summarize  the  above  analysis  by  defining  the  following  algorithm, 

<K-X0M)  Xn+1=M(^-) 

\  max  / 

where: 

Xn  =  nth  iteration  to  the  eigenvector  X 

SOLUTION  OF  THE  SYSTEM  OF  FINITE  ELEMENT  EQUATIONS 

To  reduce  storage  requirements  and  insure  computational  efficiency,  a  stand  alone  computer 
program  was  written  employing  the  present  analysis. 

The  system  of  equations  for  the  finite  element  representation  can  be  partitioned  into  a 
sequence  of  submatrices  which  may  be  expressed  in  the  following  form, 


A 

V 

B 

c 

dT 

x2 

r2 

dt 

\ 

< 

x3 

> 

“< 

R3 

N 

N 

• 

• 

s 

\ 

• 

• 

V* 

N 

\  \ 

0 

/ 

Q 

/ 

/ 

XL 

rl 

_ 

_l 

where: 

A,  B,  C  and  D  are  constant  submatrices 
R  ^ ,  R2, .  .  .  R  |_  given  column  vectors 
Xi, . . .  X|_  unknown  column  vectors 

We  shall  adopt  a  solution  procedure  which  removes  the  subarrays  which  lie  below  the  main 
diagonal.  Further,  the  technique  shall  yield  submatrices  which  are  symmetric  along  the  main 
diagonal. 

Let  us  consider  the  solution  of  the  following  set  of  equations,  i.e., 

SjXj  +  DTXj+1  -Rj 

OX,  +  CXj+1  +  DTXj+2  =  Ri+1 

where: 

Rj  =  column  vector  which  will  be  defined  later. 

Sj  =  symmetric  matrix  which  will  be  defined  later. 
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Multiplying  the  first  set  of  equations  by  DSj  ^  and  subtracting  the  result  fro  *  the  second  set, 
we  obtain  the  following  equivalent  set  of  equations, 

SjXj  +  DTXj+1  =  Rj 

0 X-,  +  (C-  DSj1  DT)  xi+1  +  DTXi+2  =  Ri+1  -  DS]1  Rj 

Since  matrices  Sj  and  C  are  symmetric  it  may  be  shown  that  the  matrix  C  -  DSj  ^  is  symmetric. 
Starting  with  the  initial  set  of  equations,  the  above  procedure  may  be  employed  to  remove  those 
arrays  which  lie  below  the  main  diagonal.  Continuing  this  process  to  the  last  set  of  equations,  we 
obtain. 


(C-  DS"^  DT)  Xl  =  Rl-  DS-^  Rl-1 

In  view  of  the  previous  work,  we  have  replaced  our  original  problem  by  that  of  determining  a  solu¬ 
tion  of  the  following  set  of  equivalent  equations,  i.e., 

AX1  +  BtX2  =  Rt 

SjXj  +  DTXj+1  =  Rj,  i  =  2 . L-1 

Slxl=RL 

where: 

S2  =  C-  BA-1  BT 
Sj  =  C  -  os:.1,  DT,  i  >  3 
Rl  =  R1 

r2=  r2-  BA"1  R-j 

Utilizing  the  method  of  back  substitution,  i.e.,  starting  with  the  last  of  the  previous  set  of  equa¬ 
tions,  we  may  write  the  solution  of  these  equations  in  the  form, 

XL  =  SL1  «L 

Xi  =  Sj‘1  (Rj  -  DTXj+1),  L  >  i  >  1 
Xi  =  A-1  (R,  -  BTX2) 

It  is  of  some  importance  to  observe  that  only  the  inverse  of  the  S  matrices  need  be  stored  in  the 
computer.  Further,  since  these  matrices  are  symmetric  only  slightly  more  than  half  the  elements 
of  these  arrays  need  be  placed  in  core.  Finally,  the  process  of  back  substitution  is  to  be  performed 
in  the  order  indicated  by  the  brackets,  i.e.,  only  vector  products  are  employed  in  the  computations. 
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REPRESENTATION  OF  THE  PHYSICAL  SYSTEM 

Since  we  have  completed  our  discussion  of  the  finite  element  model  we  may  consider  the 
problem  of  obtaining  a  mathematical  representation  of  the  physical  system. 

To  ascertain  the  motion  of  the  vibrating  system  we  shall  neglect  damping  and  assume  a  normal 
expansion,  i.e., 


w  (x,  y,  t)  =  £  aj  Wj  (x,  y)  sin  co  j  t  (8) 

i  =  1 

where: 

aj  =  magnitude  of  normalized  mode 
Wj  (x,  y)  =  i^  mode  shape. 

Differentiating  the  previous  expression,  equation  (8),  with  respect  to  t  and  evaluating  the  result  at 
t  =  0,  we  determine  that. 


I7  (x,  y,  o)  =  £  cjj  a;  Wj  (x,  y) 


at 


i  =  1 


Multiplying  this  last  result  by  Wjphdxdy  and  integrating  over  the  surface  of  the  plate,  we  obtain 
after  algebraic  reduction. 


m  —  .  .  a  w 

J7h"i,0-0|3r 


(o,  o,  o)  +y  J  Wj|S  {x,  v,  o)  dx  dy  = 


00  (  j^j _  _  Ca  f b _ 

zL  <*>  i  aj  j  4-^  wj  (o-  o)  wi  (o'  o)  +  J  Jq  wj  Wj  dx  dy 


(9) 


Recalling  that  the  initial  impact  event  is  constrained  to  a  neighborhood  of  the  origin,  we  may  write, 


fQ  Jq  wj  I?  (x'  v'0,  dxdY  =  e1  e2Wj  (0,0)  |^(o,  0,0) 


where: 


e1'  e2  =  t*ie  'en9t*1  and  width  of  the  impact  zone,  respectively. 

Since  linear  theory  implies  the  output  response  varies  linearly  with  velocity  we  consider  the  plate 
response  corresponding  to  a  unit  velocity.  Therefore,  utilizing  a  unit  velocity  in  the  conservation  of 
momentum  equation,  we  have, 

M*  1  =  |M  +  4phei  e2  j  (o,  o,  o) 


TO 
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Substituting  the  last  two  results  into  equation  (9)  we  discover  that, 


oo  ,  Cj  ra  /*b  \  jjj 

jC  wi  ai  {4^  (0'  0)  {°'  0)  +JQ  JQ  *j  dx  dY  /  =  47h  wj  (0-  o) 


Utilizing  the  orthogonality  relations  between  the  modes,  we  conclude  that. 


ai  =  ‘ 


47hwi(0'0> 


{47hwi2<0'°)  +  ^^wpdxdy. 


(10) 


Observe  that  equation  (10)  implies  that  to  ascertain  the  amplitude  of  a  specific  mode  we  must 
possess  the  value  of  the  integral  of  the  square  of  the  displacement  of  the  corresponding  mode.  To 
obtain  this  quantity,  a  numerical  approximation  to  the  integration  was  performed  in  the  computer 
program.  Such  work  is  summarized  in  appendix  C. 


Finally,  equation  (8)  was  empirically  altered  in  accord  with  the  results  of  beam  theory,  by  the 
inclusion  in  each  mode  of  an  exponential  damping  factor,  i.e., 


,  ,  _ ,  .  -0.5K  W:  t  . 

w(x,  y,  t)  =  2j  ajWj(x,  y)  e  *  sinwjt 


i=  1 


(ID 


where: 


K  =  damping  constant 

Equation  (11)  served  as  the  model  by  which  the  other  functions  of  displacement,  i.e.,  strain,  were 
similarly  modified. 

Work  performed  during  the  analysis  phase  led  to  the  creation  of  three  distinct  programs.  First, 
a  computer  program  was  developed  which  generated  initial  estimates  of  the  first  two  natural  fre¬ 
quencies  of  a  clamped  plate.  Generally,  the  approximation  procedure  yielded  first  frequencies 
which  were  approximately  ten  percent  too  large  and  second  frequencies  which  were  thirty  percent 
too  great.  Specifically,  it  was  necessary  to  multiply  the  first  approximate  frequencies  by  0.92  and 
the  second  approximate  frequencies  by  0.76  before  introducing  these  values  into  the  second  pro¬ 
gram.  These  factors  are  not  exact  and  could  change  with  the  application,  indeed,  their  only  pur¬ 
pose  is  to  improve  the  rate  of  convergence  of  the  second  program.  Finally,  we  observe  that  the  pro¬ 
gram  requires  little  CPU  time  on  the  CDC  computer  to  accomplish  its  objectives. 

Next  a  computer  program  was  developed  to  perform  the  analysis  indicated  by  the  finite  ele¬ 
ment  model.  This  program  utilized  both  the  bending  element  and  the  solution  technique  outlined 
previously  in  this  paper.  The  program  admits  a  10  by  10  node  subdivision  of  the  plate  and  since 
there  are  three  degrees-of-freedom  per  node  the  program  formulates  a  model  of  slightly  less  than 
three  hundred  equations  with  a  corresponding  set  of  unknowns.  Finally,  we  note  that  the  program 
required  a  maximum  of  slightly  more  than  four  seconds  of  CPU  time  on  the  CDC  computer  for  a 
single  problem  as  compared  to  an  approximate  thirty  three  second  CPU  time  on  the  CDC  computer 
for  a  NASTR  AN  solution  of  the  corresponding  problem. 
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The  last  computer  program  utilized  some  of  the  results  of  the  second  computer  program  to 
predict  the  maximum  strain  in  the  plate.  Values  of  damping  corresponding  to  those  adopted  in  the 
simply  supported  plate  analysis  were  employed  in  the  present  work.  Again  we  observe  that  the  pro 
gram  required  minimal  CPU  time  on  the  CDC  computer  to  perform  its  calculations. 

CONCLUSIONS 

All  three  programs  were  utilized  to  predict  the  response  of  the  clamped  plates  tested  by  Chou 
and  Flis  in  reference  (4).  No  difficulties  were  experienced  in  the  operation  of  the  programs  and  the 
results  of  this  work  are  summarized  in  figures  1  to  16  and  tables  I,  II  and  III. 

Comparison  of  the  B-series  experimental  results  and  the  theoretical  calculations  indicate  funda 
mental  differences  between  the  analysis  and  this  data.  Unfortunately,  the  experimental  results  lie 
above  the  theoretical  curve  corresponding  to  the  case  of  no  damping.  At  no  time  during  any  of  our 
previous  analysis  of  beams  and  plates  did  a  set  of  experimental  data  exceed  the  "no  damping"  com 
putations.  Thus,  it  is  difficult  to  reach  a  practical  conclusion  at  the  present  time.  Either  the  theoret¬ 
ical  model  does  not  correspond  to  the  physical  experiments,  i.e.,  rebound  of  the  impactor  before 
the  maximum  strain  occurs,  or  the  experimental  tests  are  much  more  subtle  than  one  would  nor 
mally  expect.  The  experimental  points  parallel  the  "no  damping"  computations  indicating  a  linear 
analysis  of  the  data  is  correct  and,  in  addition,  the  exx/v  and  eyy/v  data  points  from  the  same  strain 
gage  occupy  the  same  position  relative  to  the  several  damping  curves.  Thus,  there  is  an  apparent 
intrinsic  consistency  between  the  experimental  and  theoretical  data. 


In  the  case  of  the  F  and  H-series  data,  the  situation  is  considerably  different.  Most  of  the  data 
is  found  to  be  scattered  between  the  theoretical  curves  corresponding  to  damping  values  of  25  x 
10  ®  and  100  x  10  ®.  Thus,  these  two  sets  of  data  exhibit  approximately  the  same  order  of  damp 
ing  as  we  discovered  for  the  simply  supported  plates,  i.e.,  reference  (3).  Since  the  identical  plates 
were  employed  for  the  two  different  sets  of  tests  it  appears  that  the  boundary  conditions  do  not 


seriously  alter  the  value  of  the  damping  constant.  Therefore,  it  would  appear  that  the  empirical 
damping  correction  is  sufficient  to  correlate  both  the  simply  supported  and  clamped  plates  and  the 
damping  constants  are  not  greatly  affected  by  the  variation  in  the  boundary  conditions.  Finally, 


we  may  record  that  the  experimental  points  exhibit  the  theoretical  trends,  indicating  a  linear 
analysis  is  valid  and,  further,  the  exx/v  and  eyy/v  experimental  data  for  the  same  gage  occupies 
the  same  position  relative  to  the  theoretical  curves  so  that  the  data  appears  to  be  internally 
consistent. 
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Table  I  —  Physical  Configuration  of  Graphite/Epoxy  Plates  Impact  Tested  in  Reference  4 


Series  Designation 
and 

Stacking  Sequence 

Plate 

Designation 

Mass 

(Kg) 

Basic  lay-up 
(B-series) 

[(±45/02>2/±45/0/903  s 

B1 

133  x  133  x  3.6 

0.099 

B2 

133x260x3.6 

0.193 

B3 

260  x  133  x  3.6 

0.193 

B4 

260  x  260  x  3.6 

0.379 

B5 

133  x  387  x  3.6 

0.287 

B6 

387  x  133x3.6 

0.287 

B7 

387  x  260  x  3.6 

0.547 

B8 

0.547 _ 

B9 

387  x  387  x  3.6 

En 

F -series 

[  ( ±45/02) 2/ ±45/0/90]  2S 

FI 

133  x  133  x  7.6 

0.237 

F2 

133  x  260  x  7.6 

0.452 

F3 

260  x  133  x  7.6 

0.452 

F4 

260  x  260  x  7.6 

0.860 

H -series 

[+45/90/-45/+22.5/-67.5/-22.5/ 
+67.5/±45/+67.5/+22.5/-67.5/ 
-22.5/±67.5/±22.5/02/±22.5]  s 

HI 

0.213 

H2 

133  x  260  x  6.8 

0.410 

H3 

0.410 

H4 

260  x  260  x  6.8 

0.772 

3 


N  ADC-8 1250-60 


Table  II  —  Experimental  and  Theoretical  Values  of  the  Strain/Unit 
Velocity  for  B  Plate  Series 


xx/V  x  10+3  sec/M 


eyy/V  x  10  ,  sec/M 


K,  Damping 
Constant 


0.228 

0.453 

0.907 

1.73 

2.45 

2.99 

1.73 

2.45 

2.99 

1.73 

0.228 

1.73 

2.45 

2.99 

0.453 

1.73 

2.45 

2.99 

0.907 

1.73 

0.228 

1.73 

2.45 

2.99 

0.453 

1.73 

2.45 

2.99 

0.907 

1.73 

2.45 

0.228 

1.73 

2.45 

2.99 

0.453 

1.73 

2.45 

2.99 

0.907 

1.73 

1.80 

1.46 

1.24 

2.48 

2.14 

1.90 

3.47 

3.11 

2.85 

1.30 

1.21 

1.03 

1.84 

1.76 

1.58 

2.68 

2.57 

2.37 

1.09 

0.90 

0.75 

1.50 

1.28 

1.17 

2.09 

1.87 

1.76 

0.92 

0.79 

0.63 

1.26 

1.12 

0.97 

1.84 

1.63 

1.49 

m 

/25  x  10-6 

/oOx  10"6 

2.62 

2.03 

1.72 

3.57 

| 

2.99 

2.66 

4.97 

4.35 

3.99 

1.67 

1.50 

1.26 

2.33 

2.19 

1.96 

3.41 

3.21 

2.95 

1.77 

1.56 

1.33 

2.47 

2.24 

2.05 

3.48 

3.27 

3.07 

1.38 

1.14 

0.87 

1.86 

1.61 

1.35 

2.68 

2.31 

2.07 
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Table  1 1  -  Experimental  and  Theoretical  Values  of  the  Strain/Unit 
Velocity  for  B  Plate  Series  (Continued) 


xx/V  x  10+3,  sec/M 


eyy/V  x  10  ,  sec/M 


0.228 

1.73 

1.43 

1.22 

1.16 

0.99 

1.73 

1.31 

1.24 

1.07 

2.45 

1.43 

1.71 

2.99 

1.41 

1.81 

0.453 

1.73 

2.02 

1.79 

1.71 

1.53 

2.60 

1.97 

1.87 

1.66 

2.45 

2.06 

2.63 

2.99 

2.21 

2.94 

0.907 

1.73 

3.06 

2.63 

2.48 

2.29 

3.96 

2.90 

2.71 

2.50 

0.228 

1.73 

1.16 

0.79 

0.73 

0.63 

1.75 

1.56 

1.48 

1.31 

2.45 

1.10 

1.86 

2.99 

1.08 

1.84 

0.453 

1.73 

1.78 

1.22 

1.13 

1.00 

2.92 

2.44 

2.32 

2.10 

2.45 

1.63 

2.92 

2.99 

1.54 

2.72 

0.907 

1.73 

2.23 

1.67 

1.59 

1.51 

4.08 

3.29 

3.23 

3.05 

0.228  1.73  0.58  0.71 

2.45  0.57 
2.99  0.64 

0.453  1.73  1.07  1.01 
2.45  1.04 
2.99  1.02 

0.907  1.73  1.42  1.36 

2.45  1.39 


0.228  1.73  0.46  0.55 

2.45  0.45 
2.99  0.47 

0.453  1.73  0.84  1.02 

2.45  0.84 
2.99  0.90 

0.907  1.73  1.33  1.49 

2.45  1.31 


1.04  1.08 
0.98 
1.00 

1.59  1.56 

1.57 

1.57 

2.40  2.13 
2.47 


0.93  0.68 
0.96 
1.00 

1.62  1.37 

1.57 

1.57 

2.51  1.98 
2.45 
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Table  III  —  Experimental  and  Theoretical  Values  of  the  Strain/Unit 
Velocity  for  F  Plate  Series 


€  xx/V  x  10+3,  sec/M 


e  yy/V  x  1 0+^,  sec/M 


K,  Damping 
Constant 


K,  Damping 
Constant 


r 

111 

h 

fa  x  Iff6 

/ 100  x  1C6 

7 

A 

/25x  10-6 

/ 100  x  10-® 

J 

1  0.228 

1.73 

0.66 

1.35 

0.95 

0.64 

0.94 

1.95 

1.29 

0.87 

2.45 

0.68 

0.95 

2.99 

0.67 

0.98 

0.453 

1.73 

0.88 

1.94 

1.47 

1.09 

1.33 

2.77 

2.01 

1.49 

2.45 

0.91 

1.35 

2.99 

0.92 

1.39 

0.907 

1.73 

1.17 

2.67 

2.21 

1.76 

1.79 

3.77 

3.03 

2.42 

2.45 

1.22 

1.84 

0.228 

1.73 

0.56 

0.85 

0.76 

0.53 

0.72 

0.99 

0.89 

0.63 

2.45 

0.55 

0.73 

2.99 

0.53 

0.74 

0.453 

1.73 

0.82 

1.38 

1.21 

0.91 

1.11 

1.74 

1.47 

1.10 

2.45 

0.82 

1.12 

2.99 

0.82 

1.13 

0.907 

1.73 

1.21 

1.98 

1.83 

1.47 

1.72 

2.48 

2.23 

1.80 

2.45 

1.19 

1.71 

0.228 

1.73 

0.43 

0.80 

0.55 

0.40 

0.72 

1.24 

0.96 

0.72 

2.45 

0.47 

0.75 

2.99 

0.44 

0.72 

0.453 

1.73 

0.72 

1.18 

0.89 

0.69 

1.27 

1.87 

1.54 

1.20 

2.45 

0.69 

1.25 

2.99 

0.74 

1.22 

0.907 

1.73 

0.93 

1.60 

1.34 

1.13 

1.68 

2.58 

2.30 

1.94 

2.45 

0.90 

1.67 

F2 


F3 
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Table  III  —  Experimental  and  Theoretical  Values  of  the  Strain/Unit 
Velocity  for  F  Plate  Series  (Continued) 
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Table  IV  -  Experimental  and  Theoretical  Values  of  the  Strain/Unit 
Velocity  for  H  Plate  Series 


e  xx/V  x  10+3,  sec/M 


eyy/V  x  10+^,  sec/M 


K,  Damping 
Constant 


K,  Damping 
Constant 


Uj 


r 

/ 

/ 

m 

/25x  10"6 

/ 100  x  10"6 

m 

m 

/l00  x  10'6 

0.228 

1.73 

0.98 

1.59 

1.09 

0.75 

1.41 

0.99 

0.68 

2.45 

0.95 

2.99 

1.00 

Q.4S3 

1.73 

1.32 

2.32 

1.69 

1.27 

1.16 

2.07 

1.53 

1.15 

2.45 

1.34 

1.10 

2.99 

1.35 

1.10 

0.907 

1.73 

1.75 

3.19 

2.53 

2.05 

1.49 

2.85 

2.29 

1.86 

2.45 

1.80 

1.63 

0.228 

1.73 

0.76 

1.02 

0.85 

0.62 

0.61 

0.82 

0.61 

0.45 

2.45 

0.81 

0.62 

2.99 

0.82 

0.62 

0.453 

1.73 

1.39 

1.58 

1.35 

1.05 

1.07 

1.25 

1.00 

0.77 

2.45 

1.42 

1.10 

2.99 

1.46 

1.14 

0.907 

1.73 

1.95 

2.15 

2.01 

1.68 

1.50 

1.68 

1.49 

1.25 

2.45 

2.04 

1.52 

0.228 

1.73 

0.88 

0.72 

0.52 

QM 

0.90 

0.79 

0.57 

2.45 

3 

2.99 

0.68 

0.453 

1.73 

1.19 

1.43 

1.18 

0.90 

1.28 

1.44 

1.25 

0.96 

2.45 

1.22 

1.29 

2.99 

1.25 

1.32 

0.907 

1.73 

1.73 

1.92 

1.76 

1.46 

1.81 

1.97 

1.86 

1.54 

2.45 

1.76 

1.85 
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Table  IV  —  Experimental  and  Theoretical  Values  of  the  Strain/Unit 
Velocity  for  H  Plate  Series  (Continued) 


(•it  points 


Comparison  Between  Theoretical  and  Experimental  Results  for  Plate  B1 


-  igure  2  —  Comparison  Between  Theoretical  and  Experimental  Results  for  Plate  B2 


test  points 


Figure  3  -  Comparison  Between  Theoretical  and  Experimental  Results  for  Plate  B3 


Comparison  Between  Theoretical  and  Experimental  Results  for  Plate  B4 


test  points 


IMPACTOR 


igure  6  —  Comparison  Between  Theoretical  and  Experimental  Results  for  Plate  B6 


teat  points 


Comparison  Between  Theoretical  and  Experimental  Results  for  Plate  B7 


Figure  8  -  Comparison  Between  Theoretical  and  Experimental  Results  for  Plate  B8 


•  test  points 


"igure  9  -  Comparison  Between  Theoretical  and  Experimental  Results  for  Plate  FI 


©  test  points 


Experimental  Results  for  Plate  F2 


IMPACTOB 
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igure  14  -  Comparison  Between  Theoretical  and  Experimental  Results  for  Plate  H2 
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APPENDIX  A 

STIFFNESS  AND  MASS  MATRICES  FOR  A  RECTANGULAR  ELEMENT 


A  summarization  of  work  which  was  performed  in  developing  a  finite  element  representation  for 
both  the  stiffness  and  corresponding  equivalent  mass  of  an  orthotropic  plate. 
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We  assume  that  the  transverse  deflection  of  the  plate  may  be  expressed  as  a  polynomial  of  the 
following  form, 


Ay 

Ax 

■ — ► 

W  =  ao  +  a-j  f  +  a2??  +  83  F2  +  a^Jrl  +  ag??2  +  ag?3 
87  F2*?  +  a8  F^2  +  agrj3  +  a iq?3^  +  an  Tv3 


(A1) 


where: 


ag,  ai , . . . ,  a-|  t  =  unknown  constants 
Ax  =  length  of  element 
Ay  =  width  of  element 


T?  = 


Ay 


Utilizing  the  previous  assumption,  equation  A-1,  it  may  be  shown  that, 

a  =  a-|  +  2a3  f  +84  q  +  Sag?2  +  2ajTv  +  ag??2  +  Sa^gF2*?  +  a-nr?3 

(A-2) 

0  =  a2  +  34?+  2ag  r?  +  ajf2  +  2ag?Tj  +  3agr?2  +  ato?3  +  3an?^2 

where: 

a  =  Wp 

0  =  w^ 

Employing  the  expressions  for  W,  a  and  0,  equations  (A-1)  and  (A-2),  to  determine  the  values  of  the 
unknown  constants  in  terms  of  the  deflections  and  rotations  at  the  four  corner  points  of  the  ele¬ 
ment,  we  find, 

W  =  (1-3?2-  Ft?-  3T?2  +  2?3  +  3?2  t?  +  3ft?2  +  irj3  -  2?3q  -  2?T73)Wa  + 

(F-  2t2-  Tv  +  T3  +  2j2v~  ?^)“a  + 

(t?-  Tv  -  2t ?2  +  2TV2  +  r?3-  Tv3)0a  + 

(3 F2  +  2F3  -  3Fi?  -  3 Tv2  +  2?3 r?  +  2?n3)Wb  + 

{.f2  +  T3  +T2V  ~  F3^) «b  +  (Fn-  2F^2  +  TrpWb  +  (A-3) 

(-  Tv  +  3f2v  +  3?J?2  -  2?3J?  -  2?t?3)Wc  + 


A-2 
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(-  X^n  +  pfjK:  +  (-  Tv2  +  fifiWc  + 

ifn  +  3t?2  -  3?2^  -  3f  rj2  _  2^3  +  2?3rj  +  2?r?3)Wd  +  (A  3) 

(f  V  ~  If2  V  +  $3  ijloy  +  (-  V2  +  f  V2  +  t?3  -  f  7}3)0d 

Note  that  subscripts  on  the  variables  W,  a  and  0  denote  the  corner  point  at  which  the  corresponding 
quantity  is  to  be  evaluated.  Thus,  in  view  of  equation  (A-3)  we  may  express  the  displacement  W  in 
the  form, 

W  =  GA  (A-4) 


where: 


G  =  matrix  relating  transverse  displacement  W  to  the  degrees  of  freedom  of  the  nodes 
A  =  vector  of  nodal  degrees  of  freedom 


The  potential  energy  of  a  finite  element  assumes  a  form  identical  with  equation  (1),  however,  it  is 
necessary  to  limit  the  integration  to  the  domain  of  the  element  to  obtain  the  correct  result.  Con¬ 
verting  the  integral  to  dimensionless  variables  and  introducing  the  expression  for  W,  equation 
(A-1 ),  into  the  potential  energy  expression,  equation  (1),  we  find. 


df  dr;  + 


PE’?fiTfe  >».  /  /  G?f  Grr 

H/  I  Gra  GF5dfd?+'^5'  °22  / 


AxAy 


0  0 


Observe  that  subscripts  appended  to  the  matrix  G  are  intended  to  denote  differentiation  with 

1  T  A 

respect  to  the  dimensionless  variable  indicated.  Since  PE  = —  A1  K  A,  the  previous  equation  implies 
that,  2 


.1  -1 

J 

0  0 


1  1 


K  =— sr  Dn  /  /  Gl?G-r?dfdn  +-?—  H  f  f 

Ax3  11  J  J  ft  ft  AxAv  J  J 


~  ~  Ax 
d?dT?  +  Ay3  °22 


ff 

.1  J 

J 

0 


/'  / 


where: 


AxAy 


g:.  g^s  d? dr? 
rjr?  ^ 


0  0 


Ghv 


(A-5) 


K  =  stiffness  matrix 


A-3 
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Further,  we  note  that  the  maximum  kinetic  energy  of  a  finite  element  also  adopts  a  form  identical 
with  equation  (5),  however,  it  is  necessary  to  limit  the  integration  to  the  domain  of  the  element. 
Therefore,  if  we  introduce  equation  (A-4)  into  the  maximum  kinetic  energy  expression,  equation 
(5),  we  conclude  that, 


KE 


max 


JL  „,2at 


tu*  A  1  ph 


J  GTG  dA  •  A 
A 


1 

Recalling  that  KEmax  =— 


o u>2  aT 


A 

MA,  the  previous  result  implies  that 


A 

M  =  ph 


/ 


G^G  dA  =  ph  Ax  Ay 


f  GTG  dfd? 
0 


(A-6) 


where: 

A 

M  =  equivalent  mass  matrix 

Equations  (A  5)  and  (A-6)  formed  the  basis  for  the  determination  of  the  stiffness  and  equivalent 
mass  matrices.  As  a  part  of  this  task,  the  G  matrix  was  introduced  into  the  two  previously  men¬ 
tioned  expressions.  Since  the  elements  of  G  include  variables  the  integrations  indicated  in  the 
definitions  of  the  stiffness  and  mass  matrices  were  performed  by  employing  numerical  integration 
techniques. 
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APPENDIX  B 

DETERMINATION  OF  APPROXIMATIONS  FOR  THE  FIRST 
AND  SECOND  NATURAL  FREQUENCIES 


A  summarization  of  calculations  which  were  employed  in  constructing  a  computer  program  to 
obtain  approximations  for  the  first  and  second  natural  frequencies  of  a  clamped  plate. 
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v 


L 


FIRST  MODE  CALCULATIONS 

For  the  problem  under  consideration,  we  shall  assume  the  first  mode  may  be  represented  in 
the  form, 


=  f (?)  gin) 


(B-l) 


where: 


f  =  n-  r2)2d+i'0f2+a1r4) 

g  =  (1  -  r?2)2  {1  +  bQ  r?2  +  b1 q4) 

Sq,  ,  bg  and  =  arbitrary  constants 

Utilizing  equation  (B-1),  the  previously  indicated  integrals  in  equation  (6)  reduce  to  the  form. 


1 0  0 

1  -1 


If 

A  •'A 


wj2  dfd tj  =  / 

yo 

f2  df 

i 

A  1 

r 1  . 

«lf  f  df  =  I 

frr df 

1 

"ifr  df  drj  =  J 

fr  df 

92  drj 


4  drj 


r1  ry  r'  r 1 

J  J  W1rjrj  df  d77  -  J  f2  df  J  g^ 


drj 


Note  that  subscripts  appended  to  the  deflection  Wi  are  employed  to  denote  differentiation  with 
respect  to  the  dimensionless  variable  indicated. 

In  view  of  the  symmetric  nature  of  the  f  and  g  functions  we  need  only  record  the  following 
results, 

f  f2  df  *  ~  jl  +  —  <a0  +  —  a*)  +  —  (33u  +  2S-8-  +  —  a?)  | 

A  J  315  j  11  0  13  1  143  0  "0-1  17  1  J 


B-2 
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■iff  {'  '  M  +  4»  (3950  *  +  9^2'} 

64  (  2  _  U  .  1  ,_~2  106  ^  „  643  „  2i) 

TV  +T  a°  +  T®1  +  ?<3ao  +  33"  aoai  +  «gai ]f 


SECOND  MODE  CALCULATIONS 
For  the  second  mode  we  shall  assume  the  deflection  is  given  by, 

^  =  (1  -  f2)2(1  -  T?2)2  (1  -  r,  f2- r2n2)(1  +  s0f2  +  sn  r4)  * 


•(1  +  tQ  i?2  +  t1  T?4) 

where: 

r1 '  r2-  s0'  S1 '  *0  anc*  *1  =  arb'trarY  constant 
We  may  express  equation  (B-2)  in  the  following  equivalent  form, 
w2  =  g(f)  h(rj)  -  p(f)  q(tj) 

where: 

g<r>  =  (i  -  k2)2  (i  - ^  r2)  (i  +  sq  r2  +  Sl  r4) 

h(r?)  =  (1  -  i?2)2  0  -  r2V2)  (1  +  tQt?2  +  t-|  t?4) 
p«d  =r!  f2d  -  j2)2n  +  sgr2  +  s,  r4> 

q(T?)  =  r2  T?2  (1  -  T72)2  (1  +  tg  t?2  +  t1  r?4) 

Thus,  introducing  this  last  result  into  the  integrals  in  equation  (6),  we  obtain. 


w; 


dA  g2  d  xj'  h2  dr?-  2/8Pdf/  hqdtj  + J' p2df  J" q2  drj 


W 2^  dA  =  J g'2  dl;j'u'26n-2jl'p'd$j'  h'q'dr?  +j'p2  dtj'q'2  drj 


B-3 


(B-2) 
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+  2.340143E-3T,  -5.683205E-3]  $1 
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+  (1.170071E-3*r1  -  2.841602E-3)  Sq 

+  (5.683205E-3*r1  -  1.704961  E-2)  sQ 
+  8.524808E-3*r1  -  3.694083E-2  J 

(g')2df 

=  |(9.791653E-3*r^  -  2.808172E-2*r1  +  2.557442E-2)  s* 

+  [(2.8081 72E-2«r^  -  7.956487 E-2* r.,  +  7.956487E-2)  sQ 
+  2.841602E-2*r^  -  2.273282E-2*i'1  -  7.388167E-2]  s1 
+  (2.557442E-2*r^  -  7.956487E-2*r1  +  1.108225E-1)sg 
+  ( 7.956487 E-2*  -  1.477633E-1*ri)  sQ 

+  1.108225E-1*r^  +  1.219047} 

(P')2df 

=  9.791653E-3*s^  +  (2.808172E-2*sQ  +  2.841602E-2)  s1 
+  2.557442E-2*Sg  +  7.956487E-2*sQ  +  1.108225E-1 
P'  g'  df 

=  -  |(9.791653E-3*r1  -  1.404086E-2)  s* 

+  [ (2.8081 72E-2*r1  -  3.978243E-2)  sQ 
+  2.841602E-2*r1  -  1.136641E-2]  s1 

B-5 
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+  (2.557442E-2*r1  -  3.978243E-2)  s2 

+  ( 7.956487 E -2 *r1  -  7.388167E-2)  sQ 
+  1.108225E-1*r.j  | 

(g")2  dr 

=  |(1.869189*r2  -  4.134532*^  +  2.740725}  s2 

+[(4.134532*r2  -  8.848751*1^  +  5.873593)  sQ 

+  3.367299*r2  -  6.427705*r1  +  1.219047]  s1 
+  (2.740725*  r2  -  5.873593* r1  +  5.485714)  s2 

+  <5.873593*r2  -  1.219047E1*r1  +  3.657142)  sQ 

+  5.4857 1 4*  r2 -3.657142  *r1  +  1.280000El} 

(p">2dr 

-  1.8691 89*s2  +  (4.1 34532 *Sq  +  3.367299)  s1 
+  2. 740725*  Sq  +  5.873593*sq  +  5.485714 

g"  p"  dr 

=  -  I  (1.869189*r1  -  2.067266)  s2 
+  ((4.134532*r1  -4.424375)  sQ 


+  3.367299*1^  -  3.213852]  $1 
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+  {2.740725*r^  -  2.936796) 

+  (5.873593*r1  -  6.095238)  Sg 
+  5.485714*^-  1. 828571 } 

ORTHOGONALITY  RELATION 

To  reduce  the  complexity  of  the  search  procedure  in  the  case  of  the  second  mode  we  utilized 
the  orthogonality  relation.  Specifically,  this  converted  the  problem  from  that  of  determining  an 
extreme  value  to  that  of  determining  an  absolute  minimum  over  the  space  under  consideration.  The 
computer  program  which  utilized  this  relation  proved  to  be  computationally  rapid  and  convenient. 

Let  us  consider  the  following  integral, 


f  W1  w2  dA  =  JloW  9(?)  d?  y*h0( r?)  h(r?)  dr? 

A 

-  / 3q(?)  P<?)  df  J  h0(r ?)  g(r?)  dr? 

where: 

g0(f)  =  0  -  f2)2  (1  +  a0f2  +  ajJ4) 

hQ(r?)  =  (1  -  r?2)2  (1  +bQt?2  +  bi  r?^) 

It  may  be  shown  that, 

J g0gdf  =  -  ld0*rt  -dj } 

J  9qP  d£  =  d0 


where: 

dQ=  (2.841602E-3  +  1. 17007  IE -3*a0  +  5.542445E-4*a1)*s1 
+  (8.524808E-3  +  2.841602E-3*a0  +  I.^CXtflE^a^so 

+  3.694083E-2  +  8.524808E-3*3q  +  2.8416026-3*^ 

d1  =  (8.524808E-3  +  2.841602E-3*a0  +  1.170071E-3*a1)*s1 

+  (3.694083E-2  +  8.524808E-3*a0  +  2.841602E-3*a1)*s0 
+  4.063492E-1  +  3.694083E-2*ao  +  8.524808E-3*a1 
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For  the  present  circumstance  the  orthogonality  relation  may  be  written  in  the  form, 


J "  w^  w< 


_ M/4  A 

w-j  W2  dA  + — - —  =  0 


Employing  the  previous  two  sets  of  equations  recorded  in  this  section,  we  obtain, 

M/4 


d0e1  r1  +  d1  e0r2  =  d1  e1  + 


phab 


where: 


eQ,  e^  represent  factors  identical  with  dp,  d-j  except  that  the  constants  ag,  ai  are  replaced 
by  bQ,  b  i  and  Sg,  s1  by  tg,  t1 . 

Finally,  we  transform  this  last  expression  to  the  following  more  computationally  convenient  form, 

r1  =  0(d1e1  +  ^')/(doei) 

M/4 

r2  =  (1  -  ©)  •  (di  e,  +  — b  )/<d1  eQ) 

where: 


©  *  arbitrary  variable 
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APPENDIX  C 


NUMERICAL  APPROXIMATION  TO  INTEGRATION 


A  record  of  the  effort  which  was  expended  in  producing  a  technique  for  the  approximate  integra¬ 
tion  of  an  arbitrary  function.  Included  in  the  description  is  a  summary  of  those  modifications  of 
the  basic  procedure  which  are  indicated  as  a  consequence  of  the  symmetries  inherent  in  the 
problem. 


C-1 
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We  shall  develop  an  approximation  for  an  arbitrary  function  by  utilizing  the  twelve  points 
indicated  in  the  figure.  For  convenience,  we  express  the  nodal  position  in  the  x  direction  in  units 
of  n  and  in  the  y  direction  in  units  of  k.  Therefore,  let  us  attempt  to  represent  the  arbitrary  func¬ 
tion  f  over  the  x-y  plane  by  the  following  polynomial. 


£  x 

where  s  = 


It  may  be  shown  that  the  integration  of  the  above  function  over  the  center  square  in  terms  of  the 
values  of  the  function  at  the  nodal  points  is  given  by, 

h  k 

/  /  fdxdy  =  ^  hk  {  14  (f1#1 +flf0+ f0/l  +f0,0> 

~(f  2,1  +  f  1 ,2  +  f2,0  +  f0,2  +  f  1  ,-1  +  f-1 ,1  +  f-1 ,0  +  f0,-1 )  } 

where  f: :  =  value  of  the  arbitrary  function  f  at  node  i,  j. 


Finite  Element  Subdivision  of  Plate 
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Next,  let  us  consider  the  computational  difficulties  associated  with  the  subdivis  -  ,  of  the 
quarter  plate  in  accord  with  the  pattern  established  in  the  attached  figure.  Since  our  basic  Integra 
tion  pattern  utilizes  nodal  points  of  the  finite  element  it  may  De  shown  that  corner  points  of  the 
finite  element  subdivision  satisfy  the  following  constraints, 

1 .  Only  one  finite  element  corner  point  occurs  at  the  four  corner  points 
of  the  quarter  plate. 

2.  Two  finite  element  corner  points  are  found  to  occur  at  each  point  lying 
along  the  exterior  edges  of  the  quarter  plate. 

3.  Four  finite  element  corner  points  coincide  with  each  interior  point  of  the 
quarter  plate. 

If  we  extend  the  above  declarations  to  account  for  the  effects  of  the  integration  pattern  in  the 
region  under  consideration,  recalling  the  symmetry  conditions  which  are  applicable  to  our  physical 
problem,  we  conclude  that, 

1 .  The  four  corner  points  of  the  quarter  plate  occur  two  times. 

2.  Finite  element  corner  points  which  lie  along  an  edge  of  the  quarter  plate 
occur  four  times. 

3.  All  finite  element  corner  points  which  lie  in  the  interior  of  the  quarter 
plate  occur  eight  times. 

Considering  this  latter  information  it  may  be  shown  that  the  application  of  the  previously  recorded 
integration  formula  to  the  integration  of  any  function  over  the  quarter  plate  is  given  by, 

/  Afc*^  =  X  {^corner  +  ^^edge  +  ^^interiorj' 
where  f  corner  =  f  evaluated  at  corner  point 
f  edge  =  f  evaluated  at  edge  point 
f  interior  =  f  evaluated  at  interior  point 
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